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ABSTRACT 



We study the dynamics of large dust grains > 1 ^m with orbits outside of the helio- 
sphere (beyond 250 AU). Motion of the Solar System through the interstellar medium 
(ISM) at a velocity of 26 km s _1 subjects these particles to gas and Coulomb drag 
(grains are expected to be photoelectrically charged) as well as the Lorentz force and 
the electric force caused by the induction electric field. We show that to zeroth order 
the combined effect of these forces can be well described in the framework of the clas- 
sical Stark problem: particle motion in a Keplerian potential subject to an additional 
constant force. Based on this analogy, we elucidate the circumstances in which the mo- 
tion becomes unbound, and show that under local ISM conditions dust grains smaller 
than ~ 100 /im originating in the Oort Cloud (e.g. in collisions of comets) beyond 10 4 
AU are ejected from the Solar System under the action of the electric force. Orbital 
motion of larger, bound grains is described analytically using the orbit- averaged Hamil- 
tonian approach and consists of orbital plane precession at a fixed semi-major axis, 
accompanied by the periodic variations of the inclination and eccentricity (the latter 
may approach unity in some cases). A more detailed analysis of the combined effect of 
gas and Coulomb drag shows it is possible to reduce particle semi-major axes, but that 
the degree of orbital decay is limited (a factor of several at best) by passages through 
atomic and molecular clouds, which easily eject small particles. 

Subject headings: comets: general - Oort Cloud - celestial mechanics - ISM: general 



1. Introduction. 



The dynamics of dust particles in the inner Solar System has been previously addressed by 



many authors in different contexts: effects of radiation forces on circumsolar (Burns et al. 1979) 



and circumplanetary ( Horanyi et al.|1992 ) grains, electromagnetic interaction of interplanetary dust 
particles with the magnetic field of the solar wind (Landgraf 2000), grain dynamics in planetary 
magnetospheres (Horanyi 1996), and so on. 
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Recent discoveries of debris disks around young stars have brought to light additional physical 
effects related to dust dynamics such as collisions between the dust particles (Stark & Kuchner 
2009) and their resonant interaction with planetary bodies leading to asymmetries and gaps in 
debris disks (Wyatt et al. 2008). 

At the same time the dynamics of dust particles in the outer Solar System (OSS) has received 
much less attention. Here the outer Solar System means the region of space outside of the bow shock 
at 250 AU (Richardson &: Stone 2009) where unperturbed inflowing material from the interstellar 
medium (ISM; see Table [I] for a summary of ISM properties) undergoes a shock transition, and 
extending all the way to the outer edge of the Oort Cloud of comets, roughly at 5 x 10 4 AU 
(Fernandez 1999). The solar wind does not penetrate into this part of the Solar System, but the 
dust grains there move through the flow of ISM material, which perturbs their orbits. This makes 
the problem of determining the orbital evolution of grains different in the Oort Cloud as compared 
to the Kuiper Belt, since the charged component of the ISM flow does not penetrate to the Kuiper 



Belt; instead, the effects of the solar wind and planetary perturbations are important there (Moro- 



Martm & Malhotra 2002). Effects of the neutral component of the ISM flow on the Kuiper Belt 



dust particles have been considered by Scherer (2000) and Pastor et al. (2010). 

Although in this work we are concerned with the dynamics of grains in the OSS and not their 
origin, we speculate that the abundance of grains should be correlated with the spatial distribution 
of larger bodies such as comets. This is because collisions between larger bodies create a fragmenta- 
tion cascade down to smaller sizes. At the moment our knowledge of spatial distribution of comets 
heavily relies on numerical simulations of Oort Cloud formation and evolution. Such calculations 



typically produce a Cloud having both an inner and an outer edge (Kaib & Quinn 2008; Dones et 



al. 2004). The location of the outer edge in simulations is between 5 x 10 — 1 x 10 AU, which is 



only a factor of a few smaller than the typical dimension of the last closed Hill surface (Antonov 



& Latyshev 1972), beyond which objects are unbound from the Solar System by the galactic tide. 



Another result from simulations ( Kaib Sz Quinn|2008 ) is that the cometary density rises towards the 
inner edge, and it is thus likely that dust production is highest there. However, the location of the 



inner edge depends on the environment in which the Solar System formed. Kaib & Quinn (2008) 



and Brasser et al. (2006) have shown that if the Solar System formed in a cluster, the location of 



the inner edge can vary from roughly 100 AU to 3000 AU depending on the stellar density in the 
cluster. Higher stellar densities help stabilize a planetesimal kicked out by the giant planets at a 
smaller value of the semimajor axis. For this paper, we take the inner edge to be at 3000 AU, and 



results by Kaib & Quinn (2009) suggest that most of the long period comets entering the Solar 



System have initial semimajor axes at this distance. This implies that the inner edge should be no 
further than this distance, although it could be closer in. 

The goal of this work is to explore the dynamics of dust grains in the OSS by carefully analyzing 
different processes affecting their motion. Possible observational manifestations of such grains may 
provide us with information about the Oort Cloud and the collisional processes in it. Another 
reason for this study is that dust produced in the OSS may help to understand the flow of big 
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interstellar dust grains recently detected by the Galileo, Ulysses, and Cassini satellites (Grun et al. 
1994; Landgraf et al. 2000; Altobelli et al. 2003) and may contribute to the flux of micro-meteoroids 
observed at Earth ( Murray et al.|[2~004 Weryk fc Brown||2004 ) . While our work was being refereed, 



we became aware of the paper by Pastor et al. (2010) which discusses similar processes in the 
Kuiper Belt. Although their work has some similarities with our study, the methods they employed 
and some of their results are different. 

This paper is organized as follows. In ^2] we analyze the importance of different forces for the 
dynamics of dust grains. In £j3] we explore the secular effects of these forces on grain motion in 
the framework of the Stark problem. In ^4] we turn to the decay of dust particle orbits caused by 
the total drag (combined effect of gas and Coulomb drag), and in ^J5]we discuss applications of our 
results for dust evolution in the OSS. Finally, in ^6] we briefly discuss the possibility that the large 
interstellar grains observed by the satellites originate in the Oort Cloud. 



2. Forces determining grain dynamics 

If the ISM were absent in the OSS, to a zeroth order approximation, the dust particles residing 
there would move around the barycenter of the Solar System (BSS) on Keplerian orbits modified by 
radiation pressure. These orbits would slowly evolve under the action of the Poynting-Robertson 
(PR) drag, the galactic tide, and close stellar passages (Heisler & Tremaine 1986). The Keplerian 
orbital period of a body with semi-major axis a is (neglecting radiation pressure) 



T * = 27r ( G^V /2 = lo6 y r «4 /2 > W 



where a 4 = a/(10 4 AU). 



The presence of the ISM flow changes this simple picture. First, it gives rise to gas drag on the 
grains simply due to collisions of neutrals with the grain surface. Second, grains are expected to be 
photoelectrically charged to a potential of several Volts (note that Scherer (2000) and Pastor et al. 
(2010) have considered the case of neutral grains only). If the ISM has some ionized component, 



this gives rise to the Coulomb drag, which is the electric analog of dynamical friction (Binney & 
Tremaine|2008 ) and is caused by the deflection of ions around a charged grain. Third, the magnetic 
field of the ISM also interacts with charged grains. For small enough particles (or dense enough 
ISM) these forces can be stronger than the gravitational attraction to the BSS making dust grains 



unbound (see £3.3). In §£2.1 2.4 we analyze the relative effect of these and other forces on the grain 



dynamics in a variety of circumstances. 

At present, the Solar System is moving at a velocity of v w ~ 26 km s _1 through the warm 
phase of the ISM which is characterized by a gas number density n(H°) = 0.2 cm -3 , temperature 
T g = 6300 K, and ionization fraction \ ~ 0.25 (Frisch et al. 2009). The strength and orientation of 
the magnetic field carried with the wind are rather uncertain, but a typical estimate is B ~ 3 — 5 fiG 
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(Opher et al. 2009). As the Solar System moves through different phases of the ISM, the properties 
of the ISM wind may change quite dramatically compared to these numbers. Thus, we need to 
separately study the effects of the ISM flow on grains in different ISM phases. We assume for 
simplicity that the relative Solar System-ISM velocity stays equal to v w ~ 26 km s -1 at all times. 
Dust grains are assumed to be spherical and to have a density of p g = 1 g cm -3 . The grain radius 
r g is a variable parameter, but in this study we will focus on the dynamics of rather large (by 
ISM standards) particles with r g > 1 /mi, because such particles have been detected in satellite 
observations (Kruger & Grun 2009). 
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Table 1: Hydrogen number density n#, temperature T, filling factor /, ionization fraction x (Draine 
2010), and derived parameters - drag factor T (equation [9]), and quantities s e , and Si (equation 
|12| ) , which are proportional to the ion and electron Mach numbers of the grain - for various ISM 
phases. 



2.1. Electromagnetic Forces 

Solar motion relative to the ISM induces an electric field E = v„, x B/c in the Solar reference 
frame, while the magnetic field strength stays essentially the same as in the ISM frame. Letting U 
be the grain potential, the electric Fe and magnetic Fb forces on a grain are 

F E = ^v w xB, F B 3v 3 xB, (2) 

c c 

where v w and v g are the wind and grain velocities in the solar frame. The grain charge q = Ur g is 
taken to be constant, although we relax this assumption in fA) 



Grains which are large enough to be only weakly affected by the ISM wind move on (perturbed) 
Keplerian orbits at speeds which are small compared to v w : 

vk n m -1/2 , s 
— ~ 0.01a 4 , (3) 

where is the velocity of a grain moving on a circular orbit. For these grains |v 9 | = |v^| <C v w 
and \Fb\ <C |Fe| allowing us to neglect the magnetic force throughout this study when considering 
the dynamics of grains decoupled from the ISM gas flow. This assumption remains valid even for 
very small grains which get entrained in the wind provided that these particles have just been 
produced in collisions of bigger bodies (collisional debris should move with velocities < vk in the 
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Solar frame) and have not had time to get accelerated by the ISM flow to speeds comparable to 
v w . Thus, in the following we will focus only on the electric component of the electromagnetic force 



acting on grains and will mention the magnetic force only when discussing particle ejection in {3.3 



Defining 6 to be the angle between the magnetic field and the wind velocity, the strength of the 
induction electric force Fe = Ur g v w B sin 9/c relative to the gravitational force F g = GM@m g /a 2 is 

^ = 0.6 a\ Ui r~j B 5 v w>2 6 pf 1 sin 9, (4) 

where U x = U/(1Y), r 9i2 = r s /(100pm), B 5 = B/(5pG), pi = p/(lg cm" 3 ), v Wj2 6 = v w /(26 km 
s _1 ), and m g = (47r /3)p g r^ is the grain mass. Because oc r g while F g oc r?, the electric force 
becomes larger than the gravitational attraction to the BSS for particles smaller than some critical 
radius r g ^ m \ n . Grains with r g < r g ^ m \ n will be swept up in the ISM flow and ejected from the Solar 
System. We shall find in £ |3.3| that an estimate of f^min good to within about a factor of two can 
be obtained by setting Fe/F 9 = 0.25, in which case we find (neglecting the ^-dependence) 

r g ,mm = 150 pm a 4 U^ 2 B^ 2 v^ 6 p{ l/2 . (5) 

For particles smaller than r 9)m ; n , we can neglect F g compared to electromagnetic forces. Then 
in the frame of the wind, a newly created (e.g. in collisions of bigger grains) particle moves with 
speed —v w since % <C This causes gyration of the particles in the frame of the wind, while 
in the Solar System frame the particle will additionally experience an E x B drift with speed v w . 
If the angle between v w and B is not small, one expects small particles to get accelerated to a 
velocity ~ v w in the Solar System frame on a length scale 



dcoupi ~ Rl = = 1.3 x 10 8 AUr^pi v w ,26 C/f 1 Bjf 1 , (6) 

u±sr g 



where Rl is the Larmor radius of the grain. 

A notable feature of the electric force acting on grains with v g <C v w is that its magnitude and 
direction are independent of either the grain's speed or its location, provided that the magnetic field 
is homogeneous on scales comparable to the size of the Solar System. This significantly simplifies 
the analysis of the grain motion in the electric force field as we demonstrate in ^3l 



2.2. Gas Drag and Coulomb Drag 



In addition to electromagnetic forces, grains in orbit around the BSS experience gas and 
Coulomb drag due to the ISM. Since the mean free path of gas molecules and ions in the ISM is 
much larger than r g , the total drag force on a spherical grain under the assumption of sticking or 
specular reflection is given by (Baines et al. 1965 Draine &; Salpeter 1979) 



(7) 
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where 

P = nkT 

is the full ISM pressure (n is the particle number density, T is the temperature) and 
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(8) 

(9) 
(10) 
(11) 
(12) 
(13) 
(14) 



The sum in equation Q runs over all particle species i, and the first and second terms are due 
to the gas drag and Coulomb drag respectively. The term | In [l + (A/2^) 2 ] in equation ^ is a 
generalization of the usual expression for the Coulomb logarithm ln(A/zj) (Binney $z Tremaine 



2008), and Zi is the charge of species i. This is necessary because for the grain sizes we consider, 



we sometimes find A < 1. Table [T] shows the values of the dimensionless factors s and J- for ions 



and electrons in various phases of the ISM. Note that according to equation ( 10 ) even for neutral 
grains (eft = 0), the drag force does not scale quadratically with grain velocity as was assumed in 
Scherer (2000) and Pastor et al. (2010); such a scaling is only valid for s>l, which is not always 
true, see Table [T] 

We estimate the total drag force on the particles by setting v = v w , which is justified as long 
as v g <C v w (certainly true for large grains, see ^2.1): 



Fdr 



F„ 



0Aa 2 4P fr-lT 20 Pi, 



(15) 



where T20 = J-/20 and Pi = nkT/(eV cm 3 ). Unlike the electric force which scales as oc r g , the 
total drag force scales as oc r 2 , meaning that there is a critical particle radius 

?9,crit = 140 /*m Ui v w<2 6 B 5 F^qPi 1 (16) 
for which the electric force equals the total drag force. 



[2.2.1 



we abandon this simplifying assumption. 
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Additionally, if the total drag force dominates the electric force, then the minimum size of a 
bound particle obtained as before by setting Fd ra _ g / F g = 0.25 is 

r g ,mm = 160 [xma\p{ x F 2 qPi- (17) 

The coupling distance to the ISM flow d coup i = m 9 ^/(2Fd rag ), defined as the length scale over 
which the work done by the total drag force (evaluated at v = v w ) is equal to the kinetic energy of 
the grain moving at speed v w , is given by 

4oupi = 9xl0 7 AUr 9 , 2 pi vl^T^-Pr 1 - (18) 

Since the different phases of the ISM are in rough pressure equilibrium one might expect that 
-Fdrag should be comparable in all phases (one notable exception being the molecular clouds in 
which P is higher than in the average ISM owing to their self-gravity) as -Fdrag is directly related to 
gas pressure P (equation [7]). In reality we find that i*d rag varies greatly depending on the phase 
of the ISM, because the dimensionless factor T varies by orders of magnitude reflecting different 
ionization levels and Mach numbers in different phases, even though the pressure is approximately 
constant between different phases, see Table [Tj For this reason, the values of r 9)Cr it and ^ g ,min also 
vary dramatically in different environments, see Figures [T] and [7j 



2.2.1. Two dynamically distinct contributions to the total drag 

Dust grains with sizes > r g m i n feel gas and Coulomb drag only as a perturbation on top of the 
dominant gravitational force. To zeroth order, they continue moving on Keplerian orbits (albeit 
with time-varying orbital elements) and the relative speed between these grains and the ISM flow 
is v = v w + vk, where vk is the Keplerian velocity of the grain. Assuming for simplicity that the 
ISM and grain properties are constant as the grain moves on its orbit, the total drag force can be 
split as 



drag 



-F(v)- 



* drag.O 

+ AF drag + O 



VK 
Vw 



(19) 



where 



drag,0 



AF 



drag 



-F(v w )—, 

Vw 

-F(v w )— + 

V w 



dF 
dv 



+ F(v v 



(20) 
(21) 



and F(v) is the amplitude of the total drag force. The term Fd ra g,o is the zeroth order component 
which one obtains by setting the orbital velocity to zero, and AFd rag is the correction which is linear 
in vk/vw- Because Fd ra g,o is independent of and the position of the dust particle, its effect on 
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grain dynamics is analogous to that of the electric force Fe, see {2.1 We study the combined effect 
of these conservative forces on the grain motion in 

The differential drag force AFd rag (due to both Coulomb and gas drag) depends explicitly on 
vk which makes its action similar to that of a frictional force. We analyze the effect that this force 
has on the semi-major axis of dust grains in <34| 



2.3. The Galactic Tide 

The importance of the galactic tide for the dynamics of bodies in the OSS has been recognized 
for a long time. Using results of Heisler & Tremaine (1986) we find the ratio of the tidal force 
^tide = to the Solar gravitational force to be 

Ftidc - 2.6 x 10- 4 z 4 a\ ( — — ^] , (22) 



F g " " ^ 4 V 0.185 M pc-3 

where z is the vertical distance between the BSS and the object (measured in the direction per- 
pendicular to the galactic plane), = z/(10 4 AU), and po = 0.185 M pc -3 is the average stellar 



density in the galactic disk near the Sun (Bahcall 1984). 



If all other perturbations were negligible, the galactic tide would be important beyond a ~ 3000 
AU, since beyond this distance, the time for the orbital elements to cycle back to their initial values 
is shorter than the age of the Solar System ( Heisler fe Tremaine| [l986). If gas drag, Coulomb drag, 



z and scales as oc r 3 , in the particle radius, whereas the electric and drag forces are independent of z 
and scale as oc r„ and oc r 2 respectively. Thus, the tidal force dominates over the electric and drag 



and the electric force are considered, though, the picture changes. The tidal force is proportional to 

forces only for large enough r g and z, and such particles are dynamically similar to comets. Smaller 
particles, or those with smaller values of z are dynamically similar to grains, and the evolution of 
their orbital elements is governed by the electric and total drag forces. 



The particle size at which -Ftide is comparable to the electric force is 



-1/2 



r - T - n 5 cm rT 1/2 /T 1/2 U 1/2 B 1/2 v 1/2 I P ° \ (23) 

r,, cnt -U.5cmz 4 p x U x B 5 v w&y QA85 M@pc - 3 J . W 

and the particle size at which -Ftide equals the total drag force is 

r — " 8 Cm ^ * 1 Wl ( o,185M eP c-3 ) 1 ' < 24 » 

Figure [l] illustrates the cutoff between particles which are dynamically similar to grains and those 
that are dynamically similar to comets for different ISM phases. 



In this study we will not consider the effects of the galactic tide, since they have already been 
investigated by e.g. Heisler & Tremaine (1986). Instead, we will limit our discussion of the non- 



- 9 - 



gravitational forces acting on dust particles with sizes smaller than rv cr it given in equations (23) 



and (24) 



2.4. Radiation Forces 

There are several ways in which radiation affects the grain motion. Solar radiation exerts the 
radiation pressure force on dust particles 

Frad= S (g) ' (25) 

where c is the speed of light and (Q) is a frequency- averaged radiation absorption and scattering 
coefficient (Burns et al. 1979), which equals 1 in the regime of geometrical optics. Because F ra< j is 
directed radially and scales in the same way as gravity, it simply modifies F g to be F 9 (l — /3) where 
/3 = F ra _&/F g . For particles larger than 1 /xm in size (3 is small (Burns et al. 1979), and as we show 
later, particles smaller than 1 /im in size are ejected by the electric and total drag forces outside 
of the heliosphere (at a > 250 AU). For that reason we simply disregard Solar radiation pressure 
in this study (if needed the contribution of the radiation pressure can be easily accounted for by 
redefining M — > M (l — (3) in all equations). 

Poynting-Robertson drag is unimportant beyond 250 AU since the decay time for a circular 
orbit is 

16irr g pa 2 c 2 

= 2.8 x 10 4 Gyv rg^pialiQ}- 1 

Although this is shorter than the age of the Solar System for 1 /xm grains inside of 10 3 AU, 
these grains are ejected from the Solar System on significantly shorter timescales by other non- 
gravitational forces as we demonstrate in ^5] 

The anisotropy of background starlight exerts a force F an i s = A7rr 2 it(Q), where u = 10~ 12 erg 



cm - 3 ( |Draine|[2oTo| is the energy density of the background starlight and A = 0.1 (jWeingartner fc: 



Draine|[200T ) is a parameter quantifying the degree of anisotropy. We estimate that 

^ = 1.3x 10- 3 r-\ P 7 l a\ ( — )( (Q), (27) 

and using equations Q and (15) we find that F an ; s is sub-dominant compared to either the total 
drag, the electric force, or both in all ISM phases. 

Finally, there is also a non-conservative drag force that arises due to the redshifting and 
blueshifting of the background starlight as the dust grain orbits the BSS, similar to the differential 



drag described in ^2.2.1 However the magnitude of this force is ~ F an i s A _1 (vk/c) and this is 
always much smaller than |AFd rag | ~ (vk/vw) ^drag- Thus, for the problem at hand we can safely 
neglect any radiation forces on dust grains. 
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Fig. 1. — Map of regions in particle radius - orbital distance parameter space in which different 
perturbing forces dominate as shown by the different colors: yellow - electric force, blue - gas drag 
and Coulomb drag combined, red - galactic tide. Different panels correspond to different ISM 
phases as labeled on the figure. The solid line gives the particle size below which ejection occurs 
(see ^3.3) and the dashed line gives the condition a = d coup \ (equations [6] and [l8]). We assume 
that the grains have a density of 1 g cm" 3 and are at a potential of 1 Volt. The magnetic field has 
a strength of 5 /xG, and the ISM properties are the same as those listed in Table [TJ 
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The Stark Problem 



In the approximation that the ISM properties and v w are held constant and v w ^> vk, the 
electric force and the velocity-independent total drag force -Fdrag.o do not depend on the orbital 
parameters and are constant over an orbit. This reduces the problem of solving the grain dynamics 
to the classical Stark problem of motion in a Keplerian potential subject to an extra force S which 
is constant in magnitude and direction. In our case, 



S — Fe + Fdrag,0 j 

has both a component Fd ra g,o parallel to v w and a component orthogonal to v v 
general S is oriented arbitrarily with respect to v w . 



(28) 



Thus, in 



The classical Stark problem in its most general setting, including the case of S = |S| > F g , has 
been previously explored analytically. It has been shown to be separable in parabolic coordinates 



by Landau & Lifshitz (1976) and by Banks Sz Leopold (1978), who studied the ionization of highly 



excited atoms by an electric field. Kirchgraber (1970) has shown that it is possible to solve it by 
using the Kustaanheimo-Stiefel (KS) transformation (Kustaanheimo & Stiefel 1965), which maps 
out the three-dimensional Keplerian problem into the four- dimensional harmonic oscillator problem. 
Unfortunately, the results of these studies are expressed in terms of integrals of Jacobian elliptic 
functions, the special functions inverse to the elliptic integrals, and are thus difficult to analyze. 



Nevertheless, we will use the existing analyses of the Stark problem in {3.3 to obtain an accurate 
description of the conditions under which particle ejection occurs. 

In most of this study, however, we are interested in the motion of grains large enough for the 
Stark force, equation (28), to be considered a perturbation on top of F g . We can quantify this 
condition as 



Sa 2 



a = 



GM Q m 



< 1, 



(29) 



where a is the instantaneous value of the semi-major axis. In this limit, we can use standard 



methods of perturbation theory to investigate the dynamics of dust particles. Previously, Mignard 



& Henon (1984) have used osculating elements (Burns et al. 1979; Murray Sz Dermott 2001) to 



solve the problem of a constant force perturbing a test particle orbiting a central mass which is 
stationary in a rotating reference frame. This problem is relevant to the study of a particle orbiting 
a planet subjected to radiation pressure from the Sun, and it is identical to the Stark problem if 
the rotation rate of the reference frame is set to zero. 



We will approach the perturbative Stark problem (a <C 1) in two steps. First, in {3.1 



we 



consider the simplified case of the planar Stark problem in which S lies in the plane of the particle 
orbit using osculating elements. This provides us with a simple qualitative picture of the secular 
evolution under the action of Stark force. We then adopt the orbit-averaged Hamiltonian formalism 
to study the more general Stark problem in which S can have any orientation with respect to the 
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orbital plane. We also compare our results to numerical orbit integrations using an integrator 
described in Appendix [Cj 



3.1. The Planar Stark Problem: perturbative approach. 



In the case of the Stark force lying in the plane of the orbit, the motion is restricted to this 
plane, and the inclination of the orbit i = 90° does not change. We choose a coordinate frame in 
this plane so that the x-axis is aligned with S and count the longitude of pericenter to from this 



direction. The orbit averaged equations for the evolution of the osculating elements (Burns et al. 



1979 Murray & Dermott 2001 ) then become 



da 
~dt 
de 
dt 
du 
~dl 



0, 



JS 



2 GM@m g 

3 JS 



2 GM Q m g e 



sin(a;), 
cos(o;), 



(30) 
(31) 
(32) 



where J is the specific angular momentum and e is the orbital eccentricity. Both J, e, and uj vary 
in time under the action of Stark force, while a does not vary on average. 



We can integrate equations (31) and (32) directly to obtain 



ecosw 



esinw 



A. 



T s tark / 

where A is the value of ecosw at time t = 0, to is the other integration constant, and 



4vrm„ GM P 



' stark 



35 



0.67 Myra 



3/2 



(33) 
(34) 



(35) 



is the timescale for the orbital elements to return to their original values (a Stark cycle), with a 
defined in equation (29). The ratio of the Stark period to an orbital period is T stllT ^/TK = (2/3)a _1 . 



Equation (33) agrees with the results of Pastor et al. (2010), and the Stark period was previously 



obtained by Mignard & Henon (1984) 



In Figure [2] we show the evolution of the orbital elements over a timescale T star k/4, which 
corresponds to a = 0.018. The agreement between the analytical theory (equations [33] and 34 
and numerical results is good. Note that in accordance with equation (130 



the semi-major axis of 

the orbit stays constant because the constant perturbing force does no work over a closed orbit. 
However, on a time scale Tk, a still experiences small oscillations with amplitude 



Aa 



4a, 



(36) 
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Fig. 2. — Evolution of (a) semi-major axis, (b) eccentricity, (c) longitude of pericenter, and (d) 
radius over a Stark period for an orbit with T^^/Tk = 37 (a = 0.018), A = e(0) = 0.75 and 
to = ^(0) = 0. Solid lines correspond to solutions obtained from numerical integration, and dashed 
lines are the analytic solutions from equations ( 33 ) and ( 34 ) . The eccentricity tends to unity when 
the angular momentum vector changes sign, and when this happens the particle can pass very close 
to the origin. 



which are due to the work done by the Stark force in the course of orbital motion. 

The eccentricity of the orbit varies through the Stark cycle from A to 1, thus allowing very 
close approaches of particles to the Sun. As we will see in the next section, such high values of e 
are a peculiarity of the planar Stark problem, but the periodic variations of e and uj are generic. 



3.2. The General Stark Problem: perturbative approach 



Following the procedure outlined in Heisler & Tremaine (1986), we treat the general Stark 
problem of an arbitrarily oriented Stark force S by orbit-averaging the particle Hamiltonian and 
determining the integrals of motion. This procedure is different from the approach of Pastor et 
al. (2010) and Mignard & Henon (1984) who used orbit averaged equations for the evolution of 
the osculating elements (Burns et al. 1979; Murray <fe Dermott 2001) to treat the general Stark 
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problem. 

If we align the z-axis with S, the Hamiltonian per unit mass becomes 

2 r m g 

Using the relations z = rsinisin(/ + a;) and r = a(l — e 2 )/(l + e cos /), where / is the true anomaly 
(varying on an orbital timescale), we can rewrite this expression as 

H = _GMq + Sa(l-e 2 ) S m(f + Lo) S mi 
2a rn g (l + e cos /) 

When a <C 1 the orbital elements vary on a timescale T sta rk> which is much longer than Tx- 
Averaging the Hamiltonian over a closed Keplerian orbit we obtain 

ii av = ae smi sinw, (39) 

2a 2m g 

By introducing the Delaunay angle-action variables (Binney &; Tremaine 2008) 

I , L = ^GMqu (40) 
u , J = VGM a(l - e 2 ) (41) 
Q , J z = J cos i (42) 

(I is the mean anomaly, f2 is the longitude of ascending node) we rewrite the time-averaged Hamil- 
tonian as 

(GM Q ) 2 3 SL 2 I J 2 I J 2 
H - = —21^ ~ 2GM^ g V 1 " L 2 V 1 " f 2 (43) 
Because the Delaunay variables are canonically conjugate and the Hamiltonian is independent of I, 
Q, and t, we immediately see that L, J z , and must be conserved, implying that the semi-major 
axis, the z-component of the angular momentum, and total energy do not change on average. The 
total angular momentum J is not conserved, so that e, i, uj and must vary under the action of 
Stark force. 

Introducing dimensionless variables 

K = v Kz = lJ r> (44) 



we find from equation (43) that 

(A#) 2 = (l-K 2 )(l-||)sin 2 ^, (45) 

is an integral of motion which we will use from now on instead of H av . It is easy to see that 
0< AH <1-K Z . 
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Conservation of AH and K z allows us to directly relate K to to. For a given value of K z 
in K — to space, the particle is confined to move on contours of constant AH. We have plotted 
the contours for K z = 0.5 in Figure 3(a) along with the numerically integrated particle trajectory 
in these coordinates. This figure clearly shows the existence of a stationary solution when AH 
takes its maximum value of 1 — K z with both K = \/K z and to = tt/2 constant in time. The 
three-dimensional trajectory of this orbit is shown in Figure [3(b)| and it is obvious that is the 
only orbital element of the stationary orbit that varies with time (equation |51|). 



We now determine the explicit time dependence of e, i, to and f2 in the general case for arbitrary 
J z , L, and AH. The evolution of to is governed by the equation 

dto _ dH^ 3 SL K 2 Z - K A 

~db ~ dJ ~ 



sin a;, 



(46) 



2 GM Q m g K 2 ^J(l-K 2 ){K 2 -K 2 Z ) 

where H av is the Hamiltonian in equation (f43|)). Using equation ( |45| ) to express sin to and dto/dK 
as functions of K we derive the following equation for the evolution of K: 



dK 

~dt ~ 2 GM Q m g 



3 SL \{K 2 ^-K 2 )(K 2 



K 2 ■ ) 



where 



K 



f max 1 
1-min J 



1 + K 2 - {AHf 



K 



1 + K 2 - (AHf 



K 2 



1/2 



(47) 



(48) 



are the maximum and minimum possible values of K 2 for given K z and AH. Integrating equation 

K 2 ■ +K 2 K 2 

mi n ' TP a v ty 



(47) one finds 



K\t) 



+ 



K 2 ■ ft 

m ^ sin A-k 



^stark / 



(49) 



2 2 

where to is the integration constant. It immediately follows from this solution that the timescale 



for the orbit to go once around a contour in K — to space (Figure 3(a) ) is T star k/2 independent of 
K z and AH. 



Since e 2 {t) = 1 — K 2 (t) and cosi = K z /K(t) (equation 44), the solution, equation (49), 
immediately gives us the time dependence of e and i. Also, equation (45) connecting to and K 
enables us to determine explicitly the evolution of to(t). 



Finally, to determine the evolution of f2, we use the equation 
dtt dH av 3 SL y/l-K 2 K z . 



dt 



8J Z 2 GM Q m g ^K 2 - K 2 K 



smuj. 



(50) 



Equation (47) allows us to convert this expression into an equation for Q(K) with the solution 

Q (K) = arctan 




K 2 

mm 



)(^ma> 

K 2 )(K 2 m 




(51) 
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Since the argument of arctangent in the numerator of equation (J51J) vanishes at K = K mm and the 
denominator equals zero for K = K max , we have Q(K max ) — £l{K m \ n ) = tt/2. Because K goes from 
-f^min to K ma ,x on a timescale of T star k/4, £1 goes from to 2ir on a timescale of T star k. Thus, all of 
the orbital elements return to their original values after a time T star k in the general case just as in 
the planar case. The stationary solution has £1 linearly growing in time, which is most easily seen 



directly from equation (50) but can also be derived by taking the limit AH — > 1 — K z , K — > \JK, 



z 



in equation (51 ). 



These results allow us to put constraints on the possible values of eccentricity and inclination 
that an orbit attains during its Stark cycle. The minimum and maximum values of K are -fT m i n and 
-Kmax> and cosi oscillates between K z /K m&x and K z /K m \ n , while the eccentricity e varies between 
(l — K^^) 1 ^ 2 and (l — K^^) 1 ^ 2 . Since K z < K mm , the eccentricity only evolves to unity if the 
Stark vector is in the plane of the orbit (K z = 0). 



Mignard Sz Henon (1984) have previously found implicit formulas for the time dependence of 
the orbital elements. Their results also apply to a rotating reference frame, but our formulae are 
more straightforward to use in the case of the Stark problem, since they give a simple, explicit 
prescription for the evolution of the orbital elements. 

Figure |4] compares the analytic expressions for K 2 {t) = 1 — e 2 (i), i(t), co(t), and Q(t) with 
numerical integrations. Good agreement between the two approaches can be seen if one disregards 
the small-amplitude oscillations of orbital elements on a time scale Tk which are not captured by 
our orbit-averaged theory. 



We have assumed in § |3.1| and 3.2 that the Stark vector is constant, but in reality it can vary 
in both magnitude and direction due to effects such as a changing relative velocity between the 
Solar System and the ISM, a variation in the particle's charge as a function of distance from the 
Sun, small scale magnetic turbulence in the ISM, etc. We study these effects in Appendix [A] and 
conclude that most likely they add details to, but do not significantly change the general picture 
presented in §§3.1|and |3.2| 



3.3. The general Stark problem: particle ejection 



When a ~ 1 our perturbative approach used in f: 3.1 -j: 3.2 breaks down. In fact, for a > 1 one 
expects the particle to become unbound since the Stark force becomes larger than F g . We now 
perform a detailed treatment of particle ejection with the goal of determining the value of a that 
yields ejection for different initial conditions. 



Kirchgraber (1971), Dankowicz (1994), and Banks & Leopold (1978) have previously made 



attempts to study ejection in the Stark problem but their results are not easy to interpret. In 



AppendixjBJwe show that using the approach of Kirchgraber ( 1971 ) it is still possible to derive a set 
of simple analytical criteria for testing whether a particle with a given initial position r° and velocity 
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Fig. 3. — (a) Contours of constant AH (labeled by the value of AH) for K z = 0.5 with numerically 
integrated orbits over-plotted (both orbits have a = 0.02). The orbits do not exactly follow a 
contour due to our orbit- averaging of the Hamiltonian. (b) Trajectory of the stationary orbit 
(a = 0.02, K z = 0.5) for which a, e, i, and oj all remain constant, while f2 increases at a constant 
rate. The stationary orbit is marked by a cross in panel (a). 



v° is bound or unbound (equations B2 - B8 ). In finding these criteria we did not use the rather 



complicated solutions for particle motion obtained in Kirchgraber (1971), but instead resorted to 
general energetic arguments applied to the particle Hamiltonian expressed in Kustaanheimo-Stiefel 
variables (Kustaanheimo & Stiefel 1965). 

The Stark force introduces a preferred direction into the otherwise isotropic Keplerian problem. 
Thus, one expects that in general the ejection criterion would depend not only on the initial particle 
distance from the BSS but also on the direction of r° and on v . Quite interestingly, we find in 
Appendix [b] that in the case of dust grains starting at rest with respect to the BSS (v° = 0) the 
ejection criterion is independent of the direction of r° and simply states that all particles with 

a > 0.25, |v°| = 0, (52) 

are going to be ejected from the system. Figure [5] illustrates this simple ejection criterion computa- 
tionally: orbits of a number of particles with randomly chosen r° and |v°| =0 have been integrated 



for ten orbital timescales each to determine whether they are bound or not; the condition (52) 
clearly works quite well. 

To study the case of a nonzero initial velocity, we performed Monte Carlo simulations where 
we initialized particles at a given radius to have a random orientation of v° and a kinetic energy 
which is a fixed fraction e of their gravitational potential energy: (v °) /2 = eGM Q /r°. The initial 
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Fig. 4.— Plots of (a) K 2 (t), (b) fi(t), (c) i(t), and (d) u(t) for an orbit with K z = 0.6, AH = 0.34, 



and a = 0.04. Dashed lines show results from orbit-averaged theory (equations 46 - [51]), and the 
solid lines show results from numerical integrations. 



semi-major axis a can then be expressed in terms of the initial radius as a = r°/2(l — e), which 
allows the calculation of a in each case. The case of zero initial velocity corresponds to e = 0. In 
all non-zero initial velocity cases simulated here we found the general analytical ejection criteria 



(B2)-flB8b to correctly predict whether particles are bound to the Solar System. 



Figure [6] shows the effect of varying e on the value of a required for ejection. In the case of 
|v°| 7^ 0, each initial particle position is characterized by some probability of ejection since whether 
a body is ejected or not depends on the orientation of v°, and we display these probabilities with 
color maps. Most notably, Figure [6] shows that for e > there exists an outer boundary beyond 
which all particles are ejected regardless of their initial velocity, an inner boundary inside of which 
all orbits are bound for e < 1, and a region in between in which orbits can be bound or unbound 
depending on the orientation of v°. As v° increases and e becomes larger, the inner and outer 
ejection boundaries become more elongated along the direction of S, and ejection can occur for 
certain orbits even if a < 0.25, whereas some orbits can remain bound even if a > 0.25. However, 
even for e as high as 0.9 the simple ejection criterion (equation [52]) still provides a good estimate 
of the radius at which particles become unbound to within a factor of two. 

In Figure [7] we use the ejection criterion from equation (52) to determine whether particles 
are bound or unbound at different separations from the BSS in different ISM phases. It clearly 
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Fig. 5. — Illustration of the ejection criterion for particles initialized at rest, showing the outcome 
of the numerical integration of an orbit (black - ejection, red - orbits bound after an integration 
time of 10 TV). The Stark force is directed along the z-axis and because of rotational symmetry 
we only need to consider particles initialized in the x-z plane. The black circle shows the radius of 
the sphere beyond which orbits are unbound as given by the ejection criterion (equation |52|), and 
the units on the axes are scaled so that the sphere has a radius of one. The numerical integrations 
agree with equation (52), since only those particles initialized outside of a radius of one are ejected 
within 10 Tk (most unbound particles are ejected within 1 Tjc, so virtually all unbound particles 
should be ejected after 10 Tk-)- 



demonstrates that independent of the inflowing ISM phase, particles as large as r g = 150 fim should 
be ejected by the electromagnetic force at a = 10 4 AU, while at a = 10 3 AU only particles with 
< 15 /im are unbound (equation [5]). Bigger grains can also be removed in passages through the 
denser phases of the ISM. In particular, passages through molecular clouds would eject bodies as 



large as 1 m in size at 10 AU due to the total drag force, in agreement with Stern (1990). 



Very small particles for which a 3> 1 should be rapidly coupled to the ISM flow by the 
non-gravitational forces. We introduce the coupling particle size r 9iCOUp i as the size for which the 
coupling distance is equal to the initial particle semi-major axis, i.e. G?coupi(^"(/,coupi) — Particles 
with r g < r 9iCOU pi rapidly attain a velocity ~ v w in the solar frame. This makes it necessary to 
consider the Stark force as a function of velocity and take the magnetic force into consideration, 
since the relation v ~ vk <C v w is no longer valid. The dependence of r gjCOUp i on a is shown in 
Figure [T] by a dashed line. It is clear that efficient coupling to the ISM on scales of order the initial 
size of the system is possible only for particles with r g considerably smaller than the size at which 
ejection occurs. Particles with r g > r 9)COUp i passing close to the Inner Solar System should have 
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Fig. 6. — Ejection in the case of non-zero initial velocity. The Stark vector points in the +z 
direction, as indicated by the arrows, so the problem is rotationally symmetric about the z-axis, 
and we confine ourselves to the x-z plane. Each panel shows the fraction of ejected particles as a 
function of starting position for a given value of e (defined in the text). Red corresponds to all 
particles initialized at a given location being ejected and blue corresponds to all particles being 
bound after 10 Tr. The scale of the x and z axes is the same in each plot and is normalized by 
the value of the semi-major axis at which a = 0.25 in the zero velocity case (e = 0). The black 
circle shows the a = 0.25 contour for the value of e with which particles are initialized (it contracts 
inwards with increasing e because of the relation a = r°/2(l — e)). 
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Fig. 7. — Map of the ejection regimes in particle radius - orbital distance space. In the purple 



region, particles are bound regardless of the ISM phase (Class C particles, {5.3). In the adjacent 
light blue region, particles will be ejected when passing through a molecular cloud, but will be 



bound in all other phases (Class B particles, {5.2). In the green region, particles will be ejected in 



atomic clouds and molecular clouds, but will be bound otherwise (Class A particles, {5.1). In the 
yellow region particles are only bound in the coronal phase (also Class A particles, £5.1). In the 
red region, particles are unbound in all ISM phases, and the boundary of this region is given by 
equation ([5]), which is shown as a black line on the plot. 



speeds considerably different from (smaller than) v w . 



4. Orbital Decay via the Differential Drag 



We now consider the effect of the differential drag force AFd ra g on the motion of dust particles. 
AFdrag explicitly depends on the Keplerian velocity of the grain, and plugging the expression for 



AFdrag from equation (19) into the equations for the evolution of the osculating elements (Burns ct 



al. 1979 Murray Sz Dermott 2001), one obtains the following equation for the orbit-averaged rate 



of change of the semi-major axis 

/ da\ 2a 

\dt I m g v w 

where we have defined 



F(v w )(l-P) + 



dF 

d/v 



v w (3 



(53) 
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and a, b, and c are orthonormal unit vectors with a pointing towards the pericenter, c aligned 
with the particle's angular momentum vector, and b = c x a (Burns et al. 1979). Equation (53) 



agrees with the results of Pastor et al. 2010 who obtained it for the special case of a drag force that 
depends quadratically on the velocity. 

The parameter (5 has the property that < (3 < 1 for any e < 1. Although f3 varies on 
timescales of order T star k because of the secular changes in eccentricity and orientation of the orbit, 



we can average the right-hand side of equation ( 53 ) over a Stark orbit to obtain 

<«(*)) 



ooe 
2 



-it 



1 



m n v. 



g u w 



dF 

dv 



(55) 
(56) 



where angle brackets with a subscript S denote time- averaging over a Stark orbit, angle brackets 
without a subscript denote time-averaging over a Keplerian orbit, and ao is the initial value of 



the semi-major axis. The validity of averaging over a Stark orbit in equation (55) hinges upon 



7 1 S> T s tark) since this is the condition for the orbit to be only slightly modified by the differential 
drag over a timescale T sta rk- This condition is satisfied in practice, since according to equations 



(35) and (56) 



T^stark 



v,„S 



<1, 



(57) 



where S > -Fd rag (S > -Fdrag if> for example, the electric force dominates the total drag force). 
Thus, as long as vk/v w <C 1, the semi-major axis evolves exponentially with an instantaneous 
decay constant given by 7. 



An immediate consequence of equation (56) is that the semi-major axis can either grow or 



decay in time depending on the value of /3 and the sign of dF/dv. For gas drag, dF/dv > for 
any v, so 7 > 0, and the orbit always decays with time. On the other hand, for pure Coulomb 
drag dF/dv < if the Mach number of the ISM flow with respect to the sound speed of the ionized 
component s > 0.968 (equations [9j-[l4j). Thus, depending on the value of /3, the semi-major axis 
can either grow or decay with time. 

An interesting question is whether differential drag can cause the expansion of particle or- 



bits under realistic ISM conditions. According to equation (56) this amounts to a calculation of 
dF/dv\ v=v w where F is the magnitude of the force due to both gas and Coulomb drag combined 
under actual ISM conditions. The results of such a calculations are presented in Figure [8] and 
clearly show that for any grain size in any of the ISM phases dF/dv\ v=Vw > 0, meaning that a 
particle's semi-major axis can only decay (7 > 0) under the action of the differential drag force. 



5. Application to Grains in the Outer Solar System 



Having explored the forces that affect the dynamics of dust grains in the OSS, we are now in 
a position to describe their long-term evolution after they have been created in collisions of larger 
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Fig. 8. — Plot of ciln Fj din v\ v=Vw as a function of grain size for various ISM phases, where F is 
the sum of Coulomb and gas drag forces. Because dlnF/dlnv\ v=Vw > over the whole range of 
grain sizes for each phase of the ISM, the differential drag always causes the semi-major axes of 



particle orbits to decay, see equation ( 53 ) . 



objects (e.g. comets). Our discussion will predominantly focus on two important questions. First, 
what is the survival time T s of a dust particle of a given size located at a given separation from the 
BSS against ejection by the ISM flow. Second, how does the semi-major axis of a particle change 
as a result of differential drag during its lifetime in the OSS? 

In addressing these questions we must realize that through its 4.5 Gyr history, the Solar 
System has sampled different phases of the ISM, which must have subjected dust grains in the 
OSS to broadly varying ambient conditions. As a result, grains which are only weakly affected by 
the ISM flow when the Solar System passes through one phase (e.g. the coronal gas) may become 
fully coupled to this flow in another phase (e.g. a molecular cloud), see Figure [TJ Therefore there 
are several possible classes into which dust grains can be separated from the point of view of their 
long-term dynamics: 



• Class A covers particles which can survive ejection in the warm and coronal phases of the 
ISM, but which will be ejected in a passage through an atomic or molecular cloud. 

• Class B involves particles which are large enough to survive ejection in atomic clouds, but 
which would be ejected in passages through molecular clouds. 

• Class C covers particles which are large enough to survive ejection in a molecular cloud. Thus 
these particles can survive ejection in any phase of the ISM. 



The assignment of a given dust grain to a particular class depends not only on the particle size 
but also on the particle distance from the BSS. Figure [7] shows the map of different Classes in the 
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T g — a space. Note that particles satisfying the condition r g < r 9jiri i n with r ff]m i n defined in equation 
([5]) do not exhibit long-term dynamics on the OSS as they get rapidly ejected by the electric force 
even in the coronal phase of the ISM. For that reason we do not introduce a separate class for these 
small grains. 



5.1. Class A 



Stern ( 1990 ) quotes 40 Myr for the average time between passages through atomic clouds. This 
is roughly survival time, T s , of Class A particles since these are ejected in passages through atomic 
or molecular clouds, but passages through molecular clouds happen on much longer, ~ 400 Myr 



(Stern 1990), timescales. 



We now study the decay in semi-major axis caused by the differential drag in the coronal and 
warm phases of the ISM in between ejection events. We consider particles at 3000 AU with p = 1 
g cur 3 and a = 0.05, which corresponds to a grain radius r g ~ 100 fim for both the coronal and 
warm phases. The reason the grain radii are the same in different phases is because the electric 
force is the dominant perturbation in both cases, and we are assuming the particles are charged 
to the same potential of 1 V. Our choice of a is close to the threshold indicated in the ejection 



criterion (52) and is motivated by the fact that differential drag is most important for the smallest 
particles. 



Equation (55) predicts roughly exponential decay of the particle orbit on some characteristic 



time scale idecay = 7 .To evaluate upper and lower limits for tdccay, we calculate 7 1 in equation 



(56) for f3 = and f3 = 1. We can immediately rule out differential drag playing an important role 
in the coronal phase, since 2.8 Gyr < idccay < 3.0 Gyr there, which is much longer than the 40 
Myr interval between ejection events. On the other hand, a similar calculation for the warm phase 
yields 27 Myr < t^, decay < 38 Myr. Since this is smaller than T s ~ 40 Myr, the semi-major axes of 
some Class A particles can decay by a factor of a few between ejection events. 

To verify our estimate for idecay in the warm phase, we perform Monte Carlo simulations of dust 
grain orbital evolution. We initialize particles to have a semi-major axis of 3000 AU, corresponding 
to the inner edge of the Oort Cloud where most of the dust should be concentrated. Without 
good knowledge of what the distribution of the other orbital elements should be, we choose the 
eccentricities uniformly on the interval (0,0.99), and randomly orient the orbit in space. The 
perturbations included in the simulation on top of the gravitational attraction to the BSS are the 
electric force, gas drag, and Coulomb drag; other perturbations forces are unimportant for particles 
of this size at this semimajor axis. 

We simulate 10 4 orbits for 40 Myr assuming particle parameters adopted for Class A (see 
above). For each run we compute the average 7 over the length of the simulation by finding a 
best exponential fit to the dependence of the semi-major axis on time. A typical example for the 
evolution of the semi-major axis of a particle having a = 0.05, but with a starting semi-major axis 
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of 10 4 AU rather than 3000 AU is shown in Figure §1. The reason for choosing 10 4 AU for the 
starting semi-major axis is to highlight the rapid oscillations happening on an orbital timescale, 
the long oscillations occurring on a Stark period timescale, and the steady decay due to differential 
drag. The resulting distribution of decay times, idecay = l~ l -> is plotted in Figure 10 1. We calculate 
from this distribution that the mean decay time for Class A grains is (tdecay) = 30.4 Myr with a 
standard deviation of 1.1 Myr. This agrees nicely with the analytical limits given above. 



Furthermore, we see from equation ( 56 ) that the limits on £decay do not depend on the semima 



jor axis. Thus, as long as a grain is not eroded or destroyed in a collision, and the limits on idecay 
are tight, 7 will be approximately constant for a grain of a given size as it decays. Of course, smaller 



and smaller particles become bound as the semimajor axis decreases, so for a = 0.05 



10 //m) 



at 300 AU in the warm phase, the limits on idecay are a mere 2.6 Myr < idecay < 4.3 Myr, whereas 
at 10 4 AU for a = 0.05 (r g 500 /xm) they are 150 Myr < tdecay < 180 Myr, substantially longer 
than the timescale between passages through atomic clouds. 

The orbital decay of Class A particles should lead to their preferential concentration at small 
radii, which may promote fragmentation of these grains in mutual collisions and the creation of 
even smaller debris particles closer to the Sun. The implications of this effect are discussed in more 
detail in £j6j 
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Fig. 9. — (a) The decay of a particle with an initial semi-major axis of 10 4 AU and a = 0.05 
over 40 Myr in the warm phase. The rapid oscillations are on an orbital timescale and the long 
oscillations are on a Stark period timescale. The dashed line is the best fit exponential to the 
decay, (b) The decay of a Class B particle. The jumps in semi-major axis occur when the particle 
is passing through an atomic cloud. The dashed line is again the best fit exponential to the decay, 
(c) Zoomed in view of one of the passages through an atomic cloud from panel (b) . 



5.2. Class B 

Class B comprises particles which are large enough to survive ejection in atomic clouds, but 
not in molecular clouds. The time between the Solar System's passages through molecular clouds 

, so this is the survival time for Class B dust grains. 

We take the parameters of a typical Class B particle to be p = 1 g cm -3 and r g ~ 0.1 cm, 



is about 400 Myr (Sterh|p90 
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since such a particle would have a = 0.05 in an atomic cloud at a = 3000 AU. Prior to ejection in 
a molecular cloud, the particle will spend a fraction /atomic of its lifetime in atomic clouds, where 
/atomic ~ 0.01 is the filling fraction of these clouds, and a fraction 1 — /atomic of its lifetime in the 
warm and coronal phases combined, see Table [TJ For simplicity, we assume that the particle spends 
its time in either the warm phase or an atomic cloud before it is ejectecQ Then, the characteristic 
decay time is t decay = [(twism(1 - /atomic) + 7atomic/atomic)]~\ where 7wism is the value of the 
decay constant 7 in the warm phase and 7 a tomic is 7 in atomic clouds. Even though /atomic <C 1 the 
contribution to the decay from atomic clouds is still important since the total drag in them is much 
larger than in the more dilute phases of the ISM. Setting as in §5.1| (3 = and j3 = 1 in equation 



( 56 ) to get upper and lower bounds on the decay time we have 205 Myr <C ^b, decay 

< 255 Myr. This 

is shorter than the T s = 400 Myr time between passages through a molecular cloud, so the orbits 
of these particles can also decay appreciably before ejection occurs. 

We checked our analytical estimates by simulating 10 4 particle orbits for 400 Myr initialized in 



the same way as in {5.1 (with Class B particle parameters adopted here). We simulate the effect of 
passages through atomic clouds as stochastic events that happen on average every 50 Myr, so that 
there are on average 8 passages through atomic clouds per simulation, and the total time spent in 
atomic clouds per simulation is on average 4 Myr, consistent with f v = 0.01. A typical example 
of the orbital evolution of a Class B dust particle along with the best exponential fit is shown in 
Figures [9b, c. 



The distribution of decay times 7 _1 is shown in Figure 10 3. The mean decay time for this 
simulation is (idecay) = 232 Myr with a standard deviation of 24 Myr in agreement with analytical 
estimates. Note that the distributions of tdecay for Class A and Class B particles look quite different. 
This is because for Class B there is an additional source of variance that is not present for Class 
A: we model the passages through atomic clouds as stochastic events, and thus, although there 
are on average 8 cloud passages per simulation, the number and timing of passages in any given 
simulation do not have to be the same. As a result, the distribution of idecay fo r Class B particles 
is closer to Gaussian. 

We also give limits on idecay for Class B particles with starting semimajor axes at 300 AU 
and 10 4 AU. For a = 0.05 (r g = 20 /mi) in an atomic cloucj^] at 300 AU the decay time is 
2.8 Myr < t_B )deC ay < 5.2 Myr, so transport is extremely important for these particles if they are 
not destroyed on even smaller timescales (e.g. by collisions). On the other hand at 10 4 AU and 
a = 0.05 (r g = 1 cm), the decay time is 2.2 Gyr < t^^ecay < 2.6 Gyr so transport is insignificant, 



2 Assuming the particle samples the coronal phase for a fraction of time /coronal ~ 0.5 would raise fdocay by 
approximately a factor of 2 because the decay in the coronal phase is negligible. 

3 The reason that particles with a — 0.05 at 300 AU in an atomic cloud (relevant for Class B) are almost same size 
as particles with a = 0.05 at 300 AU in the warm phase (relevant for Class A) comes from the fact that the electric 
force starts to dominate the total drag in both environments for particles of size ~ 10 ^tm, and we have assumed that 
particles are charged to 1 Volt in every environment. 
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since passages through molecular clouds happen on shorter timescales. 
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Fig. 10. — (a) The distribution of decay times obtained from simulations for Class A particles 
(£5.1). Particles with p = 1 g cm -3 and r g ~ 100 /im were initialized with a semi- major axis of 



3000 AU and an eccentricity uniformly selected on the interval (0,0.99). The simulation was run 
for 40 Myr, and a decay constant was obtained by fitting an exponential to the a(t) dependence 



as shown in Figure [9l (b) The distribution of decay times for particles with p = 1 g cm and 



0.1 cm (Class B parameters, see {5.2). Eccentricities and semi-major axes are initialized in 



the same way as in (a), but the integration is now for 400 Myr and we have included the effects of 



passages through atomic clouds (see £5.2 for more details). 



5.3. Class C 

Class C particles survive passages through molecular clouds which implies that they have 
T s = 4.5 Gyr, i.e. the age of the Solar System. At 3000 AU a particle with p = 1 g cm -3 and 
r g ~ 30 cm has a = 0.05 in a molecular cloud so we take these to be the parameters of a Class C 
particle. Using values for the filling fractions from Table [T] and the formula 




where the sum runs over the different phases of the ISM, and ji and fi are the corresponding 
orbital decay constants and filling factors, the decay time is 33 Gyr < idecay < 56 Gyr, which is 
considerably longer than T s . A particle with a = 0.05 in a molecular cloud would need to have a 
semimajor axis of 1000 AU for the decay time to equal the age of the solar system (the size of such 
a particle would be r g « 3 cm). Thus, beyond 1000 AU, differential drag is unimportant for Class 
C particles, and they should remain at essentially the same semi-major axis at which they were 
created. 



5.4. Comparison with previous work 



Stern (1990) has previously estimated the perturbing force on particles in the Oort Cloud. 
However, he only considered gas dynamical drag and neglected both the Coulomb drag and the 
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electric force. Moreover, his estimates for the gas drag rely on the assumption that v w is much 
larger than the thermal velocity of the species in the wind, i.e. s»l. From Table [T] we see that 
this assumption is questionable for the warm phase and incorrect for the coronal phase. 



Stern ( 1990 ) has also considered the importance of a passing shock from a nearby supernova for 
ejecting particles and has found that a supernova at 40 pc can eject particles of radius ~ 100 /im. 
Assuming a local supernova rate of 0.02 yr _1 , the average time between supernova explosions within 
40 pc is ~ 40 Myr. This is similar to the timescale for passages through atomic clouds, and because 
a typical class A particle has a radius of r g ~ 100 //m, the effect of supernova explosions is simply 
to reduce the average lifetime of the smaller Class A particles by about a factor of two. 



Finally, Stern ( 1990 ) has pointed out that Oort Cloud objects should be eroded by high velocity 
impacts with interstellar grains. This is most relevant for Class A grains since they will be eroded 
within 40 Myr, according to Stern's estimates. Thus it may be possible that these particles are 
eroded to the point of ejection before their semi-major axes can decay appreciably. However this 
statement is speculative, because the erosion rate is not well understood and depends on both the 
composition of the grains and their 3-dimensional structure. 



6. Application to Satellite Observations 

Dust experiments on board the Ulysses and Galileo satellites have detected a significant flux of 
dust particles which appear to be co-moving with the ISM flow through the Solar System: within 
the experimental uncertainties, which are quite significant^ dust grain velocities agree both in 
direction and magnitude with v^ . 

A surprising feature of this dust component is that it contains a significant fraction of large 
grains with masses in excess of > 3 x 10~ 13 g, i.e. above the upper cutoff of the standard MRN 
dust size distribution (Mathis et al. 1977; Weingartner & Draine 2001) expected for the local ISM 



(Landgraf et al. 2000). Draine (2009a) has shown that if these massive particles are pervasive 
throughout the ISM, then (1) the mass locked up in these large grains may be inconsistent with the 
amount of refractory material potentially available for dust grain formation in the Galaxy and, (2) 
even more importantly, the reddening curve calculated accounting for the large grain population is 
grossly inconsistent with observations. This essentially excludes the idea that large grains can be 
broadly spread through the Galaxy in the amounts measured by the satellites. At the same time, 
it has proven to be difficult to identify a mechanism that would cause a concentration of such large 
grains even locally in the nearby ISM (Draine 2008). 

Another possibility for the origin of massive grains is that they come from the OSS, e.g. get 
produced in cometary collisions beyond the heliopause and become coupled to the ISM flow. This 



4 The accuracy of arrival direction determination is set by the opening angle of the dust detector, which is 140° 
for Ulysses; particle speed can be determined with a la accuracy of ~ 7 — 8 km s _1 (Grun et al. 1994). 
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possibility has been first proposed by Frisch et al. (1999), who dismissed it based on the supply 



argument: given the measured mass flux in large grains and assuming a typical radius of the dust 
production region 

(i?Oort = 5 x 10 4 AU in |Frisch et aT| fll999 |)) the whole mass of the Oort Cloud 
Mo OT t = 20 — 40 M(q would be ground down in collisions on a timescale 



, -^Oort 

t0ort ~ f ~ 10 yr ' 



(59) 



where / m = 2 x 10 20 g cm 2 s 1 is their assumed value for the flux of interstellar dust grains. 
However it is unlikely that main dust production region in the Oort Cloud lies at 5 x 10 4 AU since the 
spatial distribution of comets is expected to be concentrated towards the Cloud's inner edge, which 



may lie at 3000 AU or possibly even closer flKaib fc Quihh||2009[ |2008[ |Dones et aL||2004D . If the 
inner edge is indeed at 3000 AU, this increases the estimate of the Cloud lifetime to ioort ~ 20 Myr, 
but still does not reconcile it with the age of the Solar System. For the dust particles detected by 
satellites to have the Oort Cloud origin and for the Cloud to have a fragmentation timescale toort 
on the order of several Gyr requires the small dust particles (which get entrained in the ISM wind) 
to be mainly produced at distances of order 500 AU from the BSS. 

Our present study reveals two physical effects which may cause the concentration of interme- 
diate size grains (particles with sizes above the ejection threshold, which can produce dust grains 
with sizes below the ejection threshold in mutual collisions) towards the Inner Solar System. The 



first effect is the secular variation of eccentricity in the Stark problem (see §^3.1 3.2) which in some 



cases leads to high values of e and allows particles to venture closer to the BSS. The second effect 
is the gradual decrease in the semi-major axes of grains caused by differential drag (see Q, which 
again makes it possible for grains to come closer to the BSS. 

To investigate whether these effects can reconcile the observed flux of large ISM grains with 
their possible origin in the Oort Cloud it is necessary to nail down the location of the inner edge, 
construct a model of collisional evolution of the Cloud, investigate the transport of ejected small 
particles towards the Inner Solar System, and so on. Given the complexity of associated modeling 
we do not pursue this effort here, leaving it for future investigation. Such a study may potentially 
set useful limits on the amount of dust material that is arriving into the Inner Solar System from the 
Oort Cloud. Its results will also be relevant for interpreting observations of interstellar meteoroids 
(sizes below 100 //m) detected in radar observations (Weryk Sz Brown||2004 Murray et al.||2004). 



7. Summary 

We have investigated the dynamics of dust particles with sizes between ~ 1 fim and several 
meters in the Outer Solar System subject to the effects of the flow of interstellar gas beyond the 
heliopause. We analyzed various forces affecting the motion of small particles and showed that the 
electromagnetic force can be quite important for dust dynamics: the electric field induced in the 
Solar System frame by the magnetic field carried with the ISM flow can be the main determinant 
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of the dynamics of dust particles which are bound to the Solar System. The magnetic force is 
important for motion of small particles which get ejected from the Solar System. These statements 
are true in particular for particles smaller than 100 /im interacting with the warm phase of the 
ISM through which the Solar System is currently passing. Drag forces - Coulomb drag against the 
ionized component of the ISM flow and gas drag - are crucial for dust dynamics while the Solar 
System is passing through molecular or atomic clouds. Radiation forces are never important for 
the dynamics of dust particles bigger than a micron in size outside of the heliosphere. 

We have demonstrated that the effect of these non-gravitational forces can be well-described 
by the classical Stark problem, and have studied the ejection conditions for dust grains using the 
approach based on the Kustaanheimo-Stiefel transformation. Based on these results, we determined 
the particle sizes and semi-major axes at which they are ejected from the Solar System in different 
ISM phases. Rare passages through molecular clouds can eject particles as large as 1 m if they 
were originally located at 10 4 AU from the BSS. 

We have also explored the motion of bigger, bound grains using a perturbative approach based 
on orbit-averaging the non-Keplerian part of the Hamiltonian. This allowed us to obtain a complete 
analytical description of the system on timescales longer than an orbital period. We have shown, in 
particular, that the eccentricity of dust particles oscillates in a regular fashion reaching high levels 
in some circumstances and allowing dust grains to explore radii much smaller than their semi-major 
axes. We have also demonstrated that the component of the drag force which depends on particle 
velocity causes the decay of particle orbits, and this may be an important effect for small grains. 

Finally, we have discussed the possible relevance of these dynamical effects for the origin of 
big grains recently discovered by the Ulysses and Galileo satellites, which appear to be flowing into 
the Solar System with the ISM gas. 

We are indebted to Bruce Draine for instigating our interest in this problem, continuous en- 
couragement, numerous discussions, and advice. We thank Scott Tremaine and an anonymous 
referee for comments that helped improve the quality of the paper. This work has been financially 
supported by the Sloan Foundation and NASA grant NNX08AH87G. 
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A. Variations on the Stark Problem 

Our study until now has explicitly assumed that the Stark force is constant both in space and 
time, which in practice means that all our results are valid if S does not vary on timescales shorter 
than ~ T s tark- We now briefly discuss whether the particle semi-major axis can vary in a systematic 
fashion if we abandon the assumption of constant S. 



This preprint was prepared with the AAS IATfrjX macros v5.2. 
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A.l. Radially-Dependent Grain Charging 



We have previously assumed that a particle's charge is fixed over its orbit. However, one may 
expect the charge to vary as a function of distance from the Sun ( Kimura Mann||1998 ) because of 
increased exposure to ionizing radiation at smaller radii. Then even if the direction and magnitude 
of the induced electric field are constant, the magnitude of the electric force and the Coulomb drag 
force will still vary as functions of distance from the Sun due to the radial dependence of the grain 
charge. 

To investigate what effect a radially-dependent charge has on the orbit, we allow the grain 
charge q to vary according to a simple prescription 



(r) 



[1 + ln(l + ro/r)] , 



(Al) 



where qo and ro are constants. Thus, the charge tends to qo for r > ro and increases slowly as the 
grain approaches the Sun. We simulate the effect of grain charging numerically and choose ro to 
be the starting semi-major axis of the grain and qo such that the grain is charged to 1 V at infinity. 



As shown in Figure 11 we find that although grain charging does lead to a systematic change in the 
semi- major axis on an orbital timescale, the changes are periodic on a, tinicsccile of order T^ark 

and 

on average the semi-major axis stays constant. This is because radially-dependent grain charging 
does no work on a particle over a closed Stark orbit and cannot affect its semi-major axis in a 
systematic way. 
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Fig. 11. — Evolution of the semi-major axis with radially-dependent grain charge q{r) = 
go [1 + hi(l + ro/r)] (large amplitude curve) and with fixed grain charge qo and a = 0.02 (small 
amplitude curve). The initial semi-major axis of both grains was chosen to be 5000 AU, and the 
other orbital parameters were chosen randomly, but are the same for both simulations. 
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A.2. Time- Varying Stark Vector 

If we allow the Stark vector to vary in magnitude and in time, it will no longer be true that 
the semi-major axis or the energy of the orbit are conserved. However, if the timescale on which S 
varies is long compared to the orbital period Tk, then adiabatic invariance arguments suggest that 
the semi-major axis should not vary in a systematic fashion. Indeed, the classical Stark problem is 
separable in parabolic coordinates ( Landau Lifshitz||1976 ) and admits three adiabatic invariants 



(Banks & Leopold 1978) which are conserved when S varies slowly. The energy of a particle can 
be expressed as a function of these invariants and the value of S, meaning that as the ambient 
magnetic field and wind velocity slowly change their magnitude and direction around some mean 
values, the semi-major axis will simply oscillate around its mean value as well. 

At the opposite extreme, S may have a stochastic component which varies rapidly on timescales 
shorter than Tk, e.g. due to the small scale magnetic turbulence or random charge fluctuations on 
the grain. In this case we expect a random walk-like evolution of the semi-major axis to take place 
even if the stochastic component of S averages to zero. A similar effect of a random walk in the 
semi-major axes of comets caused by random stellar passages through the OSS has been previously 
explored in Heisler & Tremaine (1986). The resulting diffusion of the initial distribution of particle 
semi-major axes will bring some particles to smaller a and make them more gravitationally bound; 
others will move to higher a, and possibly become unbound. We expect the smallest bound dust 
grains with large initial semi-major axes to be most affected by diffusive evolution from a rapidly 
fluctuating Stark vector. 

We do not attempt to model the random walk evolution here because of the large uncertainties 
in the input ingredients that such a model will involve: spectrum and strength of magnetic turbu- 
lence, random fluctuations of the grain charge, etc. We leave this subject for future investigation. 
We do point out, though, that diffusion caused by a rapidly fluctuating Stark vector would add 



on top of diffusion caused by random stellar perturbations. Kaib & Quinn (2008) have performed 
simulations which include the effects of stellar perturbations and have found that diffusion is inef- 
fecient inside of 10 AU; only 10 — 15% of objects having a > 2 x 10 4 AU arrived there by diffusion 
from initial orbits having a < 10 4 AU after a time of 4.5 Gyr. Thus, diffusion caused by stellar 
perturbations alone will not cause significant transport in semimajor axis at the inner edge of the 
Oort Cloud, where most of the dust is likely to be concentrated. 



B. Ejection Criterion 

To understand the ejection conditions in the Stark problem, we use the analytical results of 
Kirchgraber (1971) and Dankowicz (1994) which rely on the use of KS transformation proposed by 
Kustaanheimo &; Stiefel (1965). This transformation relates six components of the 3-dimensional 
Cartesian coordinates x and velocities v to 8 components of 4-dimensional generalized coordinate 
vector u and velocity w via a non-linear prescription, see e.g. Kirchgraber (1971). Upon trans- 
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formation to KS variables, the standard Keplerian problem reduces to the problem of harmonic 
oscillator motion in four dimensions, which is easier to treat in some cases. 

In the case of the Stark problem, one can show that determining whether a given particle 
is bound or not to a gravitating center is equivalent to studying the problem (Kirchgraber 1971; 
Dankowicz 1994) of motion in i?-space (R = u\ + u\ = (r + xi)/2, u\ and U4 are the first and the 
fourth components of generalized coordinate vector u) with zero energy 



in the potential (Kirchgraber 1971) 

V{R) 



{R') 2 + V(R) = 



-2SR 3 + AhR 2 - 8ER + 4f 



(Bl) 



(B2) 







where, based on the initial conditions for x 1 
direction of the Stark vector which has components S 
for various constant factors entering this potential: 

(^2^3 V 3 X 2) — a J 



^0 ^Ch 



1) ^2) -^3/ 



and v° = (vi,V2,v®) (axis 1 is in the 
(5,0,0)) we have the following expression 
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(B3) 

(B4) 
(B5) 



Constants j and h are integrals of motion arising from the conservation of the z-component of the 
angular momentum and energy. The constant E is also an integral of motion which is generic for 



Stark problem (Landau & Lifshitz 1976). It can be rewritten as 

GM - e 



E 



1 



2 



(B6) 



x 2 + £3 is the component of r perpendicular to z-axis (or axis 1), and e z is the z- 



where p 

component of the eccentricity vector e = v x (r x v 
the pure Keplerian problem. 



GMr/r - an integral of motion unique to 



One can show that whenever the physical distance r becomes infinite (the particle is unbound) , 
the coordinate R also becomes infinite (the opposite is also true - finite R corresponds to finite r), 
which makes finiteness of R an indicator of the boundedness of particle. For this reason we only 
need to determine under which conditions there exist bound states for zero-energy motion in the 
potential V(R) — their existence and particle location in one of these states would be equivalent 
to it being bound. 



For the problem of motion in a cubic potential represented by equations (B1)-(B2) at R > 



(only positive R is physically interesting) bound states appear whenever V{R) has 3 roots all of 
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which are positive - then there is a valley in the potential curve in which the motion is bound. 
This is possible if the two extrema of V(R), a minimum and a maximum occur at non-negative R- 
and R + , and V(R-) < while V(R+) > 0. It is easy to see that 



R± 



2 h 
35 



1± 



ES 



(B7) 



To be in the valley, the particle also needs to have R° < R+. Thus, a particle is bound if and 
only if for given x , v° we have 



R- > 0, R + > 0, V(R-) < 0, V(R+) > 0, and R° = r ° + x ° < 



(B8) 



In practice, for given x , v one needs to (1) calculate the constants j,E,h using equations (B3) 



(B5), (2) calculate R- and i?+ using equation (B7), and (3) check whether all of the conditions 



(B8) are satisfied. If they are, then the motion is bound, but if at least one of them is violated 



then the motion is unbound. 



(B5) to find the ejection criterion for particles with zero 

E = 



t>j = ^2 = i>3 = 0. In this case, the constants of the motion are j = 



We can directly use equations ( B2 ) 
initial velocity 

-(r° + x°){Sr°(r - x°) - 2GM]/16r°, h = (Sz? + GM/r°)/2. Setting V(R) = to find the 
turning points of the motion and solving the resulting quadratic equation for R, the condition for 
ejection becomes (r ) 2 > GM/S. Since the initial semi- major axis a = r°/2 in the zero velocity 



case, ejection must occur for a > 0.25, as given by equation (52). 



C. Numerical Integrator 



To assist our study of orbital motion, we have written an integrator based on Burdett-Heggie 



(BH) regularization. We start from the BH regularized equation of motion (Binney &; Tremaine 



2008) 



2Er -e + r 1 



(CI) 



o?r 

ds 2 m 

where ds = dt/r, e is the eccentricity vector, E is the energy of the orbit, and F is the perturbing 



force, which need not be small. We implement equation ( CI ) in a manner similar to leapfrog ( Binney 



fc Tremaine|2008p , so our algorithm consists of a series of D 1 / 2 KD 1 / 2 steps, where K denotes a full 
kick step and D]/ 2 denotes a half drift step. Each drift step corresponds to analytically integrating 



equation (CI) without the perturbing forces for a duration of regularized time As/2, and thus is 
simply the solution to an initial value problem for a forced harmonic oscillator. Each kick step is 
then simply an instantaneous application of the perturbing forces. 



Our algorithm bears a similarity to the Wisdom-Holman integrator ( Wisdom Sz Holman|1991 ), 



and in particular to the KS regularized Wisdom-Holman integrator of Mikkola (1997). Preto & 
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Tremaine (1999) have also devised a class of algorithms with an adaptive timestep for separable 



Hamiltonians, and the integrator in this class which has a timestep proportional to 1/r, the same 
as regularization, follows a Keplerian orbit exactly except for errors in the phase. The reason for 
developing a new algorithm is so that we can easily treat the magnetic force on particles. Although 
we have found that the magnetic force is dynamically unimportant for bound particles, because it 
is much smaller than the electric force, it is important for ejected particles, and a future study of 
the trajectories of these particles would need to include this force. This rules out the algorithm of 



Preto & Tremaine ( 1999 ), since the Hamiltonian is no longer separable in the presence of a magnetic 



field. On the other hand, it should be possible to modify the algorithm of Mikkola (1997) to include 



magnetic fields, but because both time and space coordinates are changed in KS regularization, 
whereas only the time coordinate is changed in BH regularization, it is simpler to treat magnetic 
fields using BH rather than KS regularization. For this reason we have decided to write our own 



algorithm rather than modifying the one of Mikkola (1997). 



One possible second order discretization for the update of the (unregularized) velocity during 
the kick step is 



At 



E + 



V+e + V- 

2c 



x B 



m 



+ v_ 



+ v. 



(C2) 



where v_ e and v +e are the velocities right before and right after the kick step, and the magnitude 
°f Fdrag is given by equation (its direction is of course opposite to the sum of orbital and wind 
velocities). Because the kick step is assumed to happen instantaneously, r is constant during the 
kick step, and we can convert between the unregularized velocity dr/dt, and the regularized velocity 
dr/ds using the simple transformation ds = dt/r. 



One complication with implementing equation (C2) is that it is implicit since v +e appears on 



the right hand side. However, in the case of a pure electromagnetic force with no drag, equation 



(C2) can be inverted to obtain an explicit scheme for the particle advance. This is known as 



Boris's algorithm and is widely used in particle in cell simulations of plasmas (Birdsall & Langdon 
2005). The kick step of Boris's algorithm can be represented as v +e = Ej^BEj^V-e, where 



E1/2 represents half of an electric acceleration step, and B is a full magnetic rotation step. This 
algorithm is explicit, second order, time-reversible, and conserves energy. 



Rather than solving the implicit equation (C2), we extend Boris's algorithm to include the 
drag force by writing v +e = F D> i/2E 1 / 2 BE 1 /2Frj ; i/2V_ e , where F Dj i/2 is half of a drag acceleration 
step. Such an algorithm is entirely explicit, and in the case where the drag force does not depend 
on the particle velocity, the expression for the kick step above is identical to Boris's algorithm with 
the substitution gE — > q~E + Fd T a,g(v w ). This is not true in the general case for which the drag force 
can depend on the orbital velocity, so we must test our algorithm to make sure that it gives reliable 
results. 



The first test we perform is the planar Stark problem, so there is a constant perturbing force 
in the plane of the orbit. The performance of a variety of modified Wisdom- Holman integrators on 
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this problem has been studied by Rauch & Holman (1999). They found that for the energy error to 



be bounded, the pericenter passage must be resolved. Since the pericenter passage lasts for a time 

3/2 i 

t p oc r p , where r p is the radius at pericenter, regularized integrators with a timestep At oc r , 
performed better than unregularized ones, when the eccentricity approached unity. 

We compare the performance of our algorithm with that of an unregularized Wisdom-Holman 
integrator and find that it is much better at conserving energy for highly eccentric orbits (Figure 
12). The reason for this is the much finer resolution of our algorithm when passing close to the 



force center, which can be seen in Figure 13 



We also test that our algorithm can capture the effects of differential drag by making sure that 
the solutions it gives are accurate and converged. To test convergence, we study the decay rate of 
a Class A particle (£ |5.1[ ) with varying resolution. We find that the change in the decay constant 
from 50 to 1000 timesteps per orbit is less than a tenth of a percent. Because we typically use 
around 1000 timesteps per orbit, we can confidently say that we are able to resolve the decay due 
to differential drag. 

To test that the code gives not only a converged, but also an accurate result for the decay, 



we compare its results against analytical estimates for stationary orbits £3.2 In that case we can 
precisely evaluate the factor 7 in equation ( |56| ) as a function of K z . Orienting the Stark vector 
along the z-axis we have the relations (v w ■ a./v w ) 2 = (sinzsinw) 2 and (y w ■ h/v w ) 2 = (sinicosu;) 2 , 

IK Z to obtain 



which we can use together with equation ( 45 ) and K 



K z 



(C3) 



From this relation for /3 and definition (56), it is straightforward to obtain 7T star k and compare this 



with 7T star k from simulations integrated for one Stark period for pure gas drag and pure Coulomb 



drag. The results are plotted in Figure 14 and show that the simulations give not only a converged 
but also an accurate result for the differential drag. 
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Fig. 12. — (a) The fractional energy error of our BH-regularized integrator (lower curve) and 
an unregularized Wisdom-Holman integrator (upper curve) for the planar Stark problem with an 
initial eccentricity of 0.5, and a = 0.02. Both simulations were run for a total of ~ 2000 steps, (b) 
The ratio of the radius to the semi- major axis as a function of time for the same simulation shows 
that when the eccentricity becomes large and the particle passes very close to the force center, the 
error in the unregularized integrator fails to stay bounded, whereas it does stay bounded for our 
integrator. 
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Fig. 13. — Trajectory of a planar Stark orbit which has the perturbing force directed along the 
x-axis, an initial eccentricity of 0.5, and a = 0.02, for two different integrators: (a) was computed 
using our BH-regularized integrator, and (b) was computed using an unregularized Wisdom-Holman 
integrator. Both orbits were computed using ~ 2000 timesteps for the whole simulation (not per 
orbit), but our integrator has much better resolution near periapse when the eccentricity is high. 




Fig. 14. — Plot of the dimensionless decay parameter 7T star k as a function of K z for the stationary 
orbit. The curves are analytical estimates for pure gas drag (dashed) and pure Coulomb drag 
(solid) obtained using equations (53), (56), and (C3). The solid points are results obtained from 
simulations and agree with analytical estimates to within about a percent or less. 



